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Abstract 



Context. High-mass X-ray binaries are among the brightest objects of our Galaxy in the high-energy domain (0.1-100 keV). Despite 
our relatively good knowledge of their basic emission mechanisms, the complex problem of understanding their time- and energy- 
dependent X-ray emission has not been completely solved. 

Aims. In this paper, we study the energy-dependent pulse profiles of the high-mass X-ray binary pulsar 4U 0115+63 to investigate 
how they are affected by cyclotron resonant scattering. 

Methods. We analyzed archival BeppoSAX and RXTE observations performed during the giant outburst of the source that occurred 
in 1999. We exploited a cross correlation technique to compare the pulse profiles in different energy ranges and developed a relativistic 
ray-tracing model to interpret our findings. We also studied the phase dependency of the cyclotron absorption features by performing 
phase-resolved spectroscopy. 

Results. The pulse profiles of 4U 01 15+63 displayed clear "phase-lags" at energies close to those of the cyclotron absorption features 
that characterize the X-ray emission of the source. We reproduce this phenomenon qualitatively by assuming an energy-dependent 
beaming of the emission from the column surface and verify that our model is also compatible with the results of phase-resolved 
spectral analysis. 

Conclusions. We showed that cyclotron resonant scattering affects the pulse profile formation mechanisms in a complex way, which 
necessitates both improvements in the modeling and the study of other sources to be better understood. 

Key words. X-rays: binaries, pulsars: individual: 4U 0115+63. 



1. Introduction 

X-ray binary pulsars (XRBPs) are among the brightest objects 
in X-rays, and they were discovered in the ea rly 70s th anks 
to the pio neering observations car ried out by iGiacconi et al.l 
d!97lh and iTananbaum et al.l dl972l) . Many XRBPs comprise a 
young neutron star (NS), endowed with a strong magnetic field 
(B ~ 10 12 G) and a supergiant or main sequence "donor" star. 
Depending on its nature, the donor star can lose a conspicuous 
amount of matter through a strong stellar wind and/or via the 
Roche lobe overflow (see e.g jFrank et al.ll2~002l) . A large part of 
this material is focused toward the NS as a consequence of the 
strong gravitational field of this object, and then threaded by its 
intense magnetic field at several thousand kilometers from the 
stellar surface. Once funneled down to the magnetic poles of the 
NS, accretion columns may form. Here, the gravitational poten- 
tial energy of the accreting matter is first converted into kinetic 
energy and then dissipated in the form of X-rays. The luminos- 
ity generated by this process can reach values of ~10 38 erg/s 
in the energy band 0.1-100 keV (see e.g.. lPringle & R ees 1972; 
iDavidson & Ostriked[l973l To a first approximation, the mag- 
netic field lines of the NS rotate rigidly with the surface of the 
star, and thus the X-ray emission emerging from the accretion 
column is received by the observer modulated at the spin period 



of the compact object, if the magnetic and rotational axes are 
misaligned. 

When averaged over many rotational cycles, the X-ray emis- 
sion of the XRBPs folded at the spin period of the NS, i.e. the so- 
called "pulse profile", shows a remarkable stability. Significant 
variations in the shape of profiles are observed when the XRBPs 
undergo transitions between different X-ray luminosity states, 
such as reported in the cases of EXO 2030+275 and V 0332+65, 
dKlochkovet al.l l2008t iTsveankov etal]|2010l and references 
therein). They have, so far, been relatively well understood in 
terms of changes in the geometry of the intrinsic beam patterrfl 
At higher luminosities than the critical value for the characteris- 
tics o f 4U 0115+63 (L x ~ 4 x 10 37 erg s" 1 . iBasko & Sunvaevl 
1976), a radiative shock forms in a relatively extended (a few 
km) accretion column, and the radiation is emitted mainly from 
its lateral walls in the form of a "fan beam". At lower luminosi- 
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1 In this paper, the "intrinsic beam pattern" is the intensity of radi- 
ation as a function of its orientation with respect to the infinitesimal 
surface element on the column surface, measured in the frame comov- 
ing with the column. The "asymptotic beam pattern" or "beam pattern" 
is an estimate of the flux integrated along the column surface as a func- 
tion of the angle formed by the column axis with respect to the line 
of sight to the observer. It is measured at the location of the observer, 
and can be obtained from the intrinsic beam pattern after accounting for 
the general relativistic effects. Only the pulse profile can be measured 
directly from the observations. 
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ties, the radiative shock is suppressed, and matter reaches the 
base of the column with the free-fall velocity. Thus, the height 
of the column can be significantly reduced, and for very weak 
sources, it reduces to a hot-spot on the NS surface. In this case, 
the X-ray radiation escapes mostly along the axis of the column 
and the system shows a typical "pencil beam". When the X-ray 
luminosity of a source is close to the Eddington limit a com- 
bination of pencil and fan beams might apply (as proposed for 
Her X-l, Leahy 2004aB). In some cases, it has also been shown 
that the em ission can originate in hollow cones (lKrausll200U 
iLeahvl 120031, and references therein), be reprocessed onto the 
NS surface to produce a soft scattering halo around the base of 
the column, and/or be shadowed b y the upper part of the accre- 
tion stream dKraus et alJll989l 12003). 

The majority of XRBPs show significant changes in the spec- 
tral energy distribution of the X-ray emission at different pulse 
phases. This is due to a variety of factors that change during the 
rotation of the star. In particular, the cross section of the scatter- 
ing between photons and electrons trapped in a strong magnetic 
field strongly depends on the angle between the photon trajec- 
tory and the magnetic field lines, but the occultation of part of 
the column during the rotation of the NS plays also a significant 
role. 

Among the different radiation processes, the cyclotron reso- 
nant scattering is probably the most crucial for the XRBPs, as it 
gives rise to the cyclotron resonant scattering features (CRSFs; 
see, e.g., [j senberg ej^alJI 1998c lArava-Gochez & Hardin glfeOOOt 
ISchonherr et alJl2007l for recent models). These features ap- 
pear in the spectra of the XRBP sources in the form of relatively 
broad (1-10 keV) absorption lines. The CRSFs, if detected, pro- 
vide a unique tool to directly estimate the NS magnetic field 
strength, because the centroid energy of the fundamental appears 
at E cyc = 1 1.6 B 12 x (1 + z) _1 keV, and higher harmonics can be 
observed at roughly integer multiples of E cyc (Here B\2 is the 
magnetic field strength in units of 10 12 G and z is the gravita- 
tional redshift in the line-forming region. lWasserman & Shapiro] 
1 1983b . Furthermore, since the intrinsic beam pattern emerges 
from the accretion column relatively close to the surface of the 
NS, the general relativity effects related to the "compactness" of 
this object (i.e., the ratio between its mass and radius) affects the 
way this radiation is perceived by a distant observer at different 
rotational phases. 

Despite the intrinsic complexity of the problem, it was soon 
realized that the joint study of the pulsed profiles and spectral 
properties of the X-ray radiation from the XRBPs at different ro- 
tational phases would have provided an unprecedented opportu- 
nity to investigate the physics of these objects and the interaction 
between matter and radiation in the presence of strong gravita- 
tional and magnetic fields. 

Early attempts to link the spectral characteristics of the ra- 
diation emerging from an accretio n column to the pu l se profi les 
of the XRBPs were presented by Meszaros & Naeell (Il985alfbh . 
Taking some of the basic interactions between matter and ra- 
diation in the column into account, these authors were able to 
produce the first energy-dependent beam patterns that could be 
used to predict the observed properties o f pulse profiles and 
spectr a of the XRBPs. A f ew ye ars l ater. iRiffert & Meszaros! 
(119881) . iMeszaros & RifferJ (U988). and lLeahv & Lil d 19951) im- 
proved these calculations by including most of the relevant rela- 
tivistic effects: light bending, gravitational red shift, and special 
relativistic beaming near the compact object. These are essential 
for comparing the model with the real observed pulse profiles 
(IRiffert et alJll993h . 



The major problem with these models is the symmetry of 
the simulated pulse profiles, whereas the real ones are asymmet- 
ric and characterized by complex structures. The symmetry is 
mainly due to the simplifications introduced in the theoretical 
treatment of the accretion process and the magnetic field geom- 
etry. A number of attempts t o use differ ent geometries have been 
discussed bv lLeahvl (1 19901 Il99lh . and Ri ffert et alJ (1 19931) . but 
they were only partially successful. The most relevant improve- 
ment was obtained by introducing a misalignment between the 
two magnetic poles on the NS surface. 

As all the physical processes and the geometrical com- 
plications described above can h ardly be tak e n into account 
within a single theoretical model, iKraus et ail d 19951) propose 
an alternative and promising way to investigate the pulse pro- 
file formation mechanism by means of a pulse decomposi- 
tion method. From a set of observed pulse profiles in sev- 
eral energy bands and luminosity states, it is possible to de- 
termine the location of the accretion columns on the NS 
and thus decouple the geometrical effects in a nearly model- 
independent wa>0. The method was successfully applied to 
a number of XRBPs for which observatio ns with sufficientl y 
high stati stics were available: Cen X-3 (Kraus et alJ U 996). 
Her X -l dBlum&Krausll2000T). EX O 2030+275 dSasaki et all 
I2010h . A 053 5+26 dCaballero et ail l2010h . V 0332+53, and 
4U 0115+63 dSasaki et alj|201 lb . The reconstructed beam pat- 
terns show that, for EXO 2030+275 and 4U 01 15+63, the beam 
patterns could be qualitatively modeled as the combination of 
the emission from a filled column (fan and/or pencil beam de- 
pending on luminosity), plus a soft scattering halo, whereas for 
A 0535+26 and V 0332+53 the presence of a hollow column 
seems to be preferred. 

The reconstruction method presented bv lKraus et alJ (jj995) 
has so far proved to be able to unveil some of the main geomet- 
rical properties of the accretion process in the brightest XRBPs. 
However, the beam patterns have only been reconstructed for 
rather large energy intervals (several keV) and thus they cannot 
be used to investigate how effects occurring in narrow energy 
ranges, e.g., the CRSFs, can affect the beam patterns. 

In this work, we begin a study of the energy dependency 
of the pulse profile in XRBPs, at energies close to that of the 
CRSFs. As a case study, we consider the source 4U 01 15+63, so 
far the o nly XRBPs that di s played CRSFs up to the sixth har- 
monics ( Hei ndl et al.l 12004 iFerrigno et alJ l2009h . A summary 
of the known properties of this source is provided in Sect. [2] 
In Sect. [3] we describe the data analysis technique and the re- 
sults of the RXTE and BeppoSAX observations of 4U 01 15+63 
performed during the giant outburst of the source in 1999. 
We find clear evidence for phase shifts in the pulse profile of 
4U 0115+63 at energies close to those of the different harmon- 
ics of the CRSF In Sect|4] we propose a simplified model to 
interpret these shifts in terms of energy dependent variations of 
the intrinsic beam pattern (see also Appendix|A]i. We discuss our 
most relevant results and draw our conclusions in Sect. [5] 



2. 4U 0115+63 

4U 01 15+63 was di scovered in 1969 du ring one of its giant type 
II X-ray outbursts (Whi tlock et al.lll989l) . On that occasion the 
X-ray luminosity of the source rose to ~ 10 37 erg/s, a value 



2 The basic assumption is the equality of the beam patterns of the two 
accretion columns. Moreover, the decomposition is unique only if the 
angle between the rotation axis and the line of sight is known, otherwise 
this angle must be assumed. 
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that is about two orders of magnitude higher than its typical qui- 
escent level. X- ray pulsations a t a per i od of Ps _ 3.6s were 
reported first by Cominskv et al. ( 1978). Rappap ort et al.l d!978l) 
measured later the orbital parameters of the system by using SAS 
data, and estimated an orbital period of P or b = 24.3 d, an eccen- 
tricity of e = 0.34, and a value of the semi-major axis of the orbit 
of ax sin i = 140.1 It - s. Here i is the unknown inclination angle 
between the normal to the plane of the orbit and the line of sight 
to the observer. 



V 



Optical and 
635 Cas to 



IR observations permitted the object 
be identified as the companion star of 
1978T). 



4U 011 5+63 (Johns et al. 1978), located at a distance of 
7-8 kpc (Negu eruela & Okazakill200ll) . An in-depth study of the 
relation between optical and X-ray emissions from the system 
showed that the type II outbursts from 4U 011 5+63 were oc- 
currin g almost regularly every three years or so (IWhitlock et al.l 
1989), and were most likely relat ed to instabilities in the 
disk due to radiative warping ([Negueruela & Okazaki 2001; 
iNegueruela et al.ll200ll) . 

In the hard X-ray domain, the continuum energy distribu- 
tion from 4U 01 15+63 could be reasonably well modeled by us- 
ing, e .g., the phenomenological model NPEX (Nakaiima et al. 
2006). Cyclotron absorption lines were observed in the spec- 



tra of this source u p the sixth harmon i cs dWheaton et ail 
White et all 11983b [Heindl et al.l [19991 ISantangelo et al 



1979; 



1999; 



Hein dl et ail 12004; Ferrigno et all 120091) . The centroid energy 
of the fundamental line is usually at ~ll-16keV, the lowest 
among those measured from the other Galactic accreting NS 
displaying evidence of cyclotron absorption lines. The analy- 
sis of the data obtained with the GINGA satellite revealed that 
both the continuum emission from the source and the centroid 
energies of the different cyclotron lines were changi ng signif- 
icantly with the NS pulse phase dMihara et all 120041) . By us- 
ing RXTE data, Nakaii ma et al.l d2006l) show that the energy 
of the lines also depends on the source overall luminosity, and 
suggest that these variations are related to the transition from 
a radiation-dominated flow at high luminosities (with a shock 
above the stellar surface) to a matter-dominated flow at low lu- 
minosities. According to this interpretation, the authors provide 
the first constraints on the height of the accretion column where 
the cyclot ron lines are generated (~0.2-0.3 times the radius of 
the NS). iTsvgankov et al.l d2007l) . They confirm that the height 
of the accretion column decreases at lower luminosities and pro- 
vide constraints on the size of the line forming region. 

The first detailed account of the broad-band, phase-resolved 
spectral properties of the X-ray emission from 4U 0115+63 
based on a physical model has been reported by Ferrign o et al.l 
(2009). These authors apply it to one BeppoSAX data set, col- 
lected during the giant outburst o f the source in 1 999, t he self- 
consistent model developed by iBecke r & Wolff] (120071) to in- 
terpret the high-energy emission from XRBPs. They find that 
the continuum emission from the source at energies lower than 
the fundamental cyclotron line (~10 keV) is roughly constant 
through the pulse phase, and this component could be described 
by a model of thermal Comptonization of a blackbody with a 
temperature of 0.5 keV and a radius comparable to the size of the 
NS (~10km). Conversely, the spectral properties of the emission 
above ~10keV can be understood in terms of Compton repro- 
cessing of the cyclotron emission. The latter is the main cooling 
channel of the accretion column in 4U 0115+63 and has been 
observed to be more prominent at the peak of the pulse. These 
results suggested that the observed radiation from 4U 01 15+63 

during periods of intense X-ray activity (Lx ~ 10 37 erg/s) is most 



likely escaping from the lateral walls of the accretion column 
and also from a scattering halo close to the surface of the NS. 
Interestingly, appl ying the beam pattern reconstruction method 
(Krausetal .1 119951 se e also Sec t. CD on th e same data also sup- 
ports this conclusion (Sasak i et al.ll201 ll) and provides the first 
constraints on the system geometrical properties (see Sect. [4]). 

3. Observations and data analysis. 

For this work, we used the observations of 4U 0115+63 per- 
formed with Bepp oSAX (als o indicated herea fter as SAX, 
Boell a et alj|1997al) and RXTE dBradt et al.ll 1993b during the gi- 
ant outburst that occurred in 1999. A log of all the observations 
is provided in Table Q] Re sults obtained from thes e data have al- 
ready be en presented b y ISantangelo et aTl d 19991) . [Heindl et al.l 
d 19991). iHeindl et all d2004l) . iNakaiima et al.l (12006I) . and 
iFerrigno et al.l d2009l) . 

BeppoSAX data were analyzed with LHEASOFT v5.2, for 
which our research group maintained a working version of 
the Saxdas software. We considered all the data obtained with 
the high-energy nar row field instruments MECS (1.5-10.5keV , 
lBoella etalJll997bl) . HP GSPC (7-44 keV. IManzo et al.lll997l) . 
and PDS (15-100keV, iFrontera et ail 119971) . All data were 
reprocessed through the standard pipeline (saxdas v.2.1). 
Background subtraction was performed using the method of the 
earth occultation for the HPGSPC, the off-source pointing for 
the PDS, and the standard calibration files for the MECS. In the 
first observation, designated (a), we lowered the earth occulta- 
tion threshold of the HPGSPC to -1° instead of the standard 
-2°, in order to properly estimate the background. The MECS 
source extraction radius was set to 8' to optimize the S/N. 

RXTE data were analyzed wit h heasoft v6.10. W e consid- 
ered data from the PCA (2- 60 keV. iJahoda et al.ll 19961) and from 
the HEXTE (1 5-250 keV, iRothschild et al.lll998l) . Good Time 
Intervals (GTI) were created by imposing that at least 30 minutes 
passed after the satellite exited the SAA, and the elevation on 
the Earth limb was at least 10°. For the PCA, we additionally re- 
quested an electron ratio <0. 15. We modeled the background us- 
ing the estimator v3.8 for bright sources, and extracted the events 
in the mode E_12 5us_64M (time resolution 122 //S and 64 energy 
channels). The instrument specific response matrices were gen- 
erated considering the corrected energy bounds and the gain cor- 
rections performed by the EDS. We estimated the background 
for the HEXTE instrument from the available off-set pointing 
and generated the corresponding response matrix with standard 
procedures. 

3.1. Timing analysis 

The timing analysis was performed exploiting the best known 
position of 4U 0115+63 (Ra, Dec= 19.630,63.740 deg, J2000) 
in the most suitable energy ranges for the different instru- 
ments on-board BeppoSAX and RXTE (1.65-8.5keV, MECS; 
8.5-25 keV, HPGSPC; 25-100 keV, PDS; 3^5 keV, PCA; 45- 
80keV, HEXTE ). All the photon arrival times were first con- 
verted to the Solar system barycenter with the tools baryconv 
(BeppoSAX) or faxbary (RXTE), and then calculated in the 
center of mass of the binary using the available orbital ephemeris 
of 4U 0115+63 dTamuraet all 1 19921) . In each observation, the 
sourc e spin period was e stimated with a phase-shift method (see 
e.g. IFerrigno et al.l2 007). and the results are reported in TableQ] 
The relatively high uncertainty affecting the period measure- 
ments is mainly due to the timing noise related to the variation 
of the pulse profile from one pulse to the other. This produces 
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Table 1. Log of the BeppoSAX and RXTE observations of 4U 0115+63 used for this work. The uncertainty on spin period 
measurements is at lcr c.l. on the last digit. The source flux is measured from the spectral analysis and the uncertainty is only 
statistical at 90% c.l. Assuming a source distance of 7 kpc, a flux of 2 x 10~ 8 erg s _1 cirT 2 would correspond to L2-iookeV = 1.2 X 
10 38 ergs -1 . 



Obs 


Instrument 


Start Time [UT] 


Stop Time [UT] 




Exposures [ks] 


P[s] 


Flux (2-1 00 keV) 




SAX 






MECS 


HPGSPC 


PDS 




[10~ 8 ergcirT 2 s~'] 




RXTE 






PCA 


HEXTE-A 


HEXTE-B 






(a) 


SAX 


1999-03-06 14:48 


1999-03-07 14:23 


54.9 


48.8 


44.3 


3.61456[2] 


1.61+0.01 


(b) 


RXTE 


1999-03-07 04:09 


1999-03-07 09:26 


17.1 


7.3 


7.2 


3.614593[6] 


1.70 + 0.02 


(c) 


RXTE 


1999-03-12 00:02 


1999-03-12 07:30 


13.8 


7.7 


7.7 


3.61428[2] 


2.05 ± 0.01 


(d) 


SAX 


1999-03-19 17:05 


1999-03-20 08:42 


31.2 


32.0 


30.0 


3.613816[3]* 


1.62 + 0.01 


(c) 


RXTE 


1999-03-21 06:27 


1999-03-21 08:33 


8.8 


1.5 


1.4 


3.61441[8] 


2.05 ± 0.04 


(f) 


SAX 


1999-03-26 17:31 


1999-03-27 17:34 


53.7 


42.5 


48.3 


3.614246[2] 


1.07 + 0.02 


(g) 


RXTE 


1999-03-28 01:04 


1999-03-28 08:12 


8.8 


1.7 


1.7 


3.61443[1] 


1.13 + 0.08 


(h) 


RXTE 


1999-03-29 04:15 


1999-03-29 11:17 


5.9 


3.5 


3.5 


3.61444[1] 


1.04 + 0.06 


* The period is referred to MJD 51 256.7037037 (TBD), when P = -(1.56 


± 0.06) x 10 


Vs. 







a scatter of ~2% on the phase of the first and second Fourier 
harmonics used in this analysis. 

To extract the energy dependent pulse profiles and the phase 
resolved spectra, we binned the events recorded in each energy 
channel of each instruments in 100 (SAX) or 128 (RXTE) phase 
bins to form the "phase-energy matrices". These matrices were 
then rebinned to achieve a sufficient S/N in each energy bir0 
(at least 180 for BeppoSAX, and 60 for RXTE ). The exposure 
time of each phase bin was computed using the GTI of the cor- 
responding observation and the reference time to fold the light 
curves was chosen so that the maximum of the high-energy pulse 
profiles in different observations remains at phase ~ 0.3. 

3.2. Determination of phase lags. 

For all the observation in Table [TJ we extracted a pulse profile 
in each of the energy bins of the phase-energy matrices (see 
Sect. 13. 11 1. The pulse profiles showed a remarkable dependence 
on the energy (especially below ~5 keV), but did not change dra- 
matically between the different observations. As an example, we 
show in Figure[TJthe pulse profiles extracted for the observation 
(a). Above ~5 keV, a prominent structure appears at phase ~0.3 
(the "main peak") and remains virtually stable at all the higher 
energies. The predominance of this structure is also confirmed 
by the analysis carried out in the left hand panel of Fig. [3] Here 
we represent with red (blue) colors the phases of all the pulse 
profiles extracted from the observations in Table [TJ where the 
source count rate was higher (lower). 

The relatively broad red spot, corresponding to the mean 
peak, seems to move in phase with the energy, thus showing 
a "phase-lag" in the pulse profiles. We note that this behavior 
is present in all observations. To investigate the origin of these 
phase lags in more detail, we extracted a reference pulse profile 
in the 8-10. 5keV energy range for each observation (see Fig. [2J 
and then computed the cross correlation between this reference 
profile and the pulse profiles extracted in different energy bands. 
We limit this analysis to the phases corresponding to the main 
peak (see Fig. fJJ and use the discrete cross correlation function: 

XiPr(i)p e (i + j) 

Cf 0) = ; „ , „ ■ (i) 

3 We optimized these values for our analysis on a trial basis. The 
different values come from the different energy resolutions of the in- 
struments. 



Here, C e (j) is the linear correlation coefficient calculated at the 
energy bin e and phase bin /'. The indices of the phase bin for 
each pulse profile are j and i, while p r and p e are the reference 
and energy dependent pulse profiles, from which we subtracted 
their average values. The maximum value (and relative phase) 
of the correlation parameter for each energy bin, C e , max , was es- 
timated as a function of the phase by performing parabolic fits. 
Hereafter, we use normalized phase units, in which the phase is 
a real number in the range (0, 1). A value of C e , max << 1 in- 
dicates that the pulse profiles p e extracted in the energy bin e 
are not strongly correlated with the reference pulse profile p r , so 
their shapes differ significantly. A value of C e max ~ 1 indicates 
that one profile can be obtained from the other using a linear 
functiorQ and intermediate values indicate a moderately differ- 
ent shape. The results of this analysis are shown in the right hand 
panel of Fig.[3]in a color representation and they provide a clear 
confirmation of significant phase lags in the energy-dependent 
pulse profiles. 

A more quantitative representation of the measured phase 
lags is provided in Fig. [4] Here, we report the phase lags of 
the main peak for each observation with respect to the refer- 
ence pulse, together with the corresponding uncertainties and the 
maximum correlation coefficient. To determine the phase lags 
and the associated uncertainties in the most appropriate way, we 
simulated 100 "phase-energy" matrices in which the number of 
events in each phase and energy bin was randomly selected from 
a Poisson distribution with an average value equal to the mea- 
sured number of counts. We then carried out a correlation analy- 
sis for each of the simulated matrices and recorded the estimated 
phase lags. The measurement reported for each observation and 
energy bin in Fig. [4] was thus determined as the average of the 
simulated values and the corresponding uncertainty taken as one 
standard deviation form this value. We verified a posteriori that 
the values measured from the real data are consistent with the 
value obtained from the average of the simulations within the 
reported uncertainties. The results reported in Fig. [4] show that 
the pattern of the phase lags vs. energy is very similar for all 
the observations. The main peak of the pulse is shifted toward 
lower phases at energies corresponding to the different absorp- 
tion cyclotron features of 4U 0115+63. In correspondence of 
each of these shifts in phase, we also notice a deviation of the 



4 The profiles in Eq. ((TJ have zero average but not unitary standard 
deviation, therefore in the ideal case of C f max = 1 , they are related by a 
function that produces a phase shift and a rescale. 
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(d) SAX 1999-03-19 17:05 - 1999-03-20 08:42 
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Figure 3. Left panel: pulse profiles extracted from all the observations in Table [T] as a function of energy. Red (blue) colors indicate 
higher (lower) count rates at different pulse phases and energies. As discussed in Sect. 13.11 the higher intensity is around phase ~0.3, 
in correspondence with the "main peak". The pulse profiles in this figure were all scaled to have in each energy bin a zero average 
value and a unitary standard deviation. Right panel: values of the correlation parameter C e , as a function of phase lag with respect 
to the reference pulse profile and energy. Red (blue) colors indicate values of C e close to unity (<sl). See text for details. 



cross correlation coefficient from unity. This testifies that, at en- 
ergies close to those of the different cyclotron features, the shape 
of the main peak is slightly different from that of the reference 
pulse profile. At the highest energies ( ~ 40 keV), the marked 
decrease in the correlation coefficient is due to statistical noise. 
To correctly evidence the features visible that are in the pattern 
of the phase lags as a function of energy, we fit a model com- 
prising several asymmetric Gaussian lines to these datfl 

We note that the intrinsic energy resolutions of the different 
instruments (not accounted for in our phenomenological fit of 
the phase lags) affects the estimate of the centroid energies of 



5 The empirical function used to fit the phase-lag energy dependency 
is /(e) = A + Eil Nj exp [- (e - E t f /crf ldM ], where a mM = a i4 for 
e < Ef, o",,(rf,„| = cr, „ for e > E h and e is the energy. The fourth Gaussian 
is only poorly determined due to the low S/N of the data. 



the Gaussian empirical functions. From Fig.|4j it is clear that the 
third Gaussian is systematically detected with a centroid energy 
lower in the PDS than in the PCA data, because the former in- 
strument is characterized by a lower energy resolution than the 
latter. The opposite effect is visible on the second Gaussian by 
comparing the PCA with the HPGSPC data, as the latter has a 
finer energy resolution than the former. Taking these systematic 
effects into account, we conclude that the phase lags measured 
from the different instruments and observations in TableQ]are in 
good agreement with each other, and no prominent variations of 
the phase lag pattern can be appreciated in the luminosity range 
of 4U 01 15+63 spanned by the data. 

Finally, we checked that all the above results do not depend 
on the particular settings adopted in the method developed for 
the analysis described above. For this purpose, we performed a 
further analysis of the observation (b) by adopting different val- 
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Figure 4. Phase lags of the main peak in normalized phase units (points) and linear correlation coefficients (stars) as a function of 
energy. We considered here BeppoSAX data from the MECS (2-8.5 keV), the HPGSPC (8.5-25 keV), and the PDS (>25keV). 
RXTE data are from the PCA (3-45 keV) and the HEXTE-B (> 45keV). The correlation function is computed on a phase interval 
of 0.4 phase units, centered on the maximum of the reference pulse profiles. The solid line represents the best-fit function to the 
phase lags (see text for the details). 



ues of the range of both the phases used for the correlation anal- 
ysis (A(p) and the energy range chosen to extract the reference 
pulse profile (AE). All the results of this analysis are reported 
in Fig. [5] We note that the adoption of a smaller A0 resulted 
in a significant decrease in the reliability of our analysis, as a 
large part of the peak falls outside the sampled window (also in- 
dicated by the lower value of the cross-correlation coefficient). 
Choosing a larger A</> would instead lead to a loss of coherence 
in the cross-correlation due to the complex energy dependency 
of the pulse profiles outside the main peak. Changing the value 
of AE resulted, as expected, only in a systematic offset of the 
phase shifts. We thus concluded that the general trend of the 
phase shifts is not affected by the particular set of parameters 
used for the analysis. 



3.3. Spectral analysis 

In order to deeply investigate the relation between the phase lags 
and the energy of the cyclotron scattering features, we also per- 
formed a spectral analysis of the source. 



We first extracted the phase-averaged spectr£@for each ob- 
servation in Table Q] and performed a fit to these data with 
the phenomenologi cal multicompone nt model proposed for 
4U 0115+63 bv iFerrigno et alJ d2009). For PCA, we added a 
0.5% systematic error and fixed the absorption column at A^h = 
9 x 10 21 crrT 2 , the best-fit value obtained in the BeppoSAX ob- 
servations. This provided, in each case, an accurate measure- 
ments of the flux (and thus the luminosity) of the source (see 
Table [TJ. For each observation, we also report the centroid en- 
ergy of the first three absorption lines in Fig. [6] together with 
the energy range at which the most negative phase lag occurs 
in the corresponding energy range. Although the most negative 
phase shifts occur slightly above the measured centroid energies 
for all lines, there is a clear indication of the link between timing 
and spectral features. In the following paragraphs, we clarify the 
origin of this apparent mismatch of energies. 

A phase-resolved spectral analysis was carried out to study 
the variation in the properties of the cyclotron absorption lines at 



For the phase-averaged spectral analysis of the RXTE observations, 
we used the standard2 mode of the PCA data and both HEXTE units. 
For the spectral fitting, we used the XSPEC software (v. 12.6.0u) 
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Figure 1. Pulse profiles of 4U 01 15+63 during observation (a) in 
different energy ranges. Data were collected with the MECS (2- 
8.5 keV), the HPGSPC (8.5-24 keV) and the PDS (24-90 keV). 
The phase-energy matrices of the three instruments were arbi- 
trarily shifted by a few energy bins (maximum 3) to take their 
timing systematics into account and achieve an optimal align- 
ment in their common energy ranges (8-1 1 keV for the MECS 
and HPGSPC, 20-25 keV for the HPGSPC and the PDS ). 



different phases (see e.g. iMihara et alj "2004). Following the dis- 
cussion in Sect. 13.11 only the range of phases including the main 
peak was considered. We extracted source spectra in 64 and 50 
phase bins for the PCA and the HPGSPC, respectively. We re- 
stricted the energy range for the phase-resolved analysis to the 
10-35 keV interval and adopted in the fit a simplified spectral 
model comprising a power-law modified by a high-energy cut- 
off and three Gaussian absorption lines (gabs in XSPEC ). This 
permitted us to clearly identify the centroid energy of the funda- 
mental and of the first harmonic in most spectra. We did check a 
posteriori that the simplified spectral model did not change our 
re sults significant l y with respect to those already reported, e.g., 
bv lFerrigno et al.l d2009l) . 

In Fig. [7] we show the values of the centroid, amplitude (cr), 
and optical depth (t) of the fundamental vs. the pulse phase, 
while in Fig. [8] we show the corresponding results for the first 
harmonic for the three BeppoSAX observations. The pulse pro- 
files in the relative energy ranges are also shown. For the sec- 
ond harmonic, we verified that the poorer S/N prevents such de- 
tailed analysis, by also using the PDS data to cover the range 



2 0.10 - 




Figure 2. Reference pulse profile of 4U 0115+63 extracted 
from observation (b). Data are from the PCA (energy range 8- 
10.5 keV). The pulse profile has been linearly scaled to have zero 
average and unitary standard deviation. The two vertical solid 
lines enclose the phase interval used for the correlation analysis: 
0.33-0.72 (the phase is shifted with respect to Fig.[3]for clarity's 
sake). 



above 34keV. The parameters of the lines could not be well con- 
strained at different phases from those of the main peak, as in 
these cases the lines appeared absent or shallower and less pro- 
nounced. Despite the slightly different shape of the pulse profiles 
and the variation in the luminosity of the source (see Table [TJ, 
the three observations gave fairly similar results for the line pa- 
rameters. This confirmed that the scattering region is located at 
roughly the same height above the NS in this luminosity range. 
In the upper right hand panel of these figures, we use a horizontal 
band to indicate the energy range in which the largest negative 
shift in phase of the pulse profiles in the BeppoSAX observa- 
tions was measured. In both cases, the centroid energies corre- 
spond to the energy of the most negative phase shift in some 
phase interval. For the fundamental, such a phase range corre- 
sponds to the highest optical depth of the line, but to an inter- 
mediate value for the first harmonic. We notice that the phase 
resolved spectroscopy proves how the phase-averaged spectral 
results are a non-trivial average of phase variability and cannot 
be directly compared to our results for phase lags. For all the 
lines, the centroid energy estimated from the average spectra is 
lower than the energy of the corresponding phase lags, but the re- 
sults of the phase-resolved spectroscopy clearly indicate that the 
energy dependency of the shift in phase of the pulse profiles is 
related to the cyclotron-scattering process. Similar results were 
obtained by performing a phase-resolved analysis of the RXTE 
data over a slightly wider range of luminosity. 

It is particularly interesting to note that in all cases the cy- 
clotron lines get narrower and deeper in the descending edge of 
the pulse (i.e., the values of cr and t are respectively lower and 
higher in these phases). We discuss this observational property 
of the line at the end of Sect. [4] 



4. A model for the emission from an accretion 
column 

In this section, we describe the model we developed to inter- 
pret the energy-dependent phase shifts of the pulse profiles from 
4U 0115+63 and the properties of its cyclotron emission lines 
(see Sect. l3.2l and l3.3l ) in terms of an energy-dependent beamed 
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Figure 5. Phase shifts measured from the observation (b) by ut 
ing a different range of phases than adopted for the analysis i 
Fig. [4] to compute the shifts (A</>) and a different energy rang 
of the reference pulse profile {AE). Specific values of these pa 
rameters are indicated in different panels. The shaded interval 
represent the centroid energy of the asymmetric Gaussians use 
to describe the lags with the l <x uncertainty. The solid line is th 
best-fit function of the phase lags. 
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Figure 6. Phase-averaged luminosity and centroid energy of the 
absorption lines. The horizontal cyan bands are the energy range 
of the most negative phase shift of the main peak measured from 
the BeppoSAX data. Uncertainties on fitted parameters are at 
90% c.l. Upper left: X-ray luminosity (see TableQ]). Upper right: 
centroid energy of the fundamental. Lower left: centroid energy 
of the first harmonic. Lower rig/jfxentroid energy of the second 
harmonic. 
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emission emerging from two accretion columns on the NS. The 
model has a two-step structure: in the first part, we compute the 
flux emerging from the column as a function of its inclination an- 
gle with the line of sight (the "asymptotic beam pattern"). In the 
second part , we exploit t h e geo metric parameters of 4U 01 15+63 
derived by ISasaki et alJ d201 ll) using the pulse decomposition 
method to model the observed pulse profiles. 

4.1. Asymptotic beam pattern of an accretion column. 

We considered in the model a spherical NS, with rad i us R^ s 
and mass Mns. Based on the conclusion of Sasa ki et aD (I201 lh . 
we exclude the "hollow cone" geometry and then assume that 
the accretion column is homogeneously filled by the accreting 
material and characterized by a conical shape. The latter is a 
good approximation in case the accr etion flow is confined by 
a dipolar magnetic field (Leahv 2003). Finally, we assume that 
the radiation is emitted from the lateral surface of the column, 
as expected for the X-ray luminosity of 4U 01 15+63 in the ob- 



Figure7. Phase-resolved study of the fundamental absorption 
line for the BeppoSAX observations. Blue points correspond 
to observation (a), green circles to observation (d), and red di- 
amonds to observation (f). Upper left: Pulse profiles extracted in 
the energy range 1 1.2-13. 2keV (shifted in phase to be aligned). 
Upper right: Centroid energy of the first cyclotron harmonic. 
The horizontal cyan band is the energy range of the most nega- 
tive phase shift of the main peak measured from the SAX data. 
Lower left: Amplitude of the line. Lower right: Optical depth of 
the line. 



Bask o & Sunvaevlll976l) . As the emission from the accretion 
column is relatively close (at most a few km) to the NS sur- 
face, we also take the most relevant relativistic effects into ac- 
count. We include in our model the relativistic beaming, the 
gravitational redshift and the gravitational lensing effect (also 
known as solid angle am plification), according to the method de- 
scribed by |Leahvl d2b03). However, for the light bend ing, we use 
servations analyzed here (Lx — 10 38 ergs _1 > 4 x 10 37 ergs~' the approximate geodesic equation proposed by iBeloborodovl 
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Figure 8. Phase-resolved study of the first harmonic for the 
BeppoSAX observations. Blue points correspond to observa- 
tion (a), green circles to observation (d), and red diamonds to 
observation (f). Uncertainties on fitted parameters are at 90% 
c.l. Upper left: pulse profiles extracted by using data from the 
HPGSPC in the energy range 21.7-23.5keV. Upper right; cen- 
troid energy of the first cyclotron harmonics. The horizontal cyan 
band is the energy range of the most negative phase shift of the 
main peak measured from the BeppoSAX data. Lower left: am- 
plitude of the line. Lower right: Optical depth of the line. 



(2002fl The mathematical details of the models are reported 
in Appendix [A] The model parameters are thus the NS mass 
and radius, the semi-aperture of the accretion column, ao, the 
column height with respect to the NS surface, and the emis- 
sivity of the surface element of the column in its rest frame, 
^o.v^, 0o, 0i) (i-e., the intrinsic beam pattern). Here, v is the 
radiation frequency, r the height at which the photon is emit- 
ted, 0o its azimuth around the column axis, a the angle between 
the local radial direction and the photon trajectory, and <p\ is the 
corresponding azimuth (see Fig. IA.ll ). 

Deriving Io tY (r,(f>o,a,(pi) on a physical basis would require 
a detailed treatment of all the relevant photon-scattering pro- 
cesses in the accretion column, which is beyond the scope of the 
present paper. Instead, we consider a phenomenological model 
of Io, v (f, (/>o,a,<pi) depending just on the angle a. The depen- 
dency on <p\ can be neglected as the cyclotron scattering in the 
geometry considered here is almost independent from this an- 
gle, because the magnetic field lines are nearly perpendicular to 
the NS surface at the polar caps. We also assume that the inten- 
sity of the radiation emitted from the column does not depend 
on the radius r and the angle (po. Including the dependence on 
0o could be relevant for tracing the coupling between the accre- 
tion disk with the magnetosphere. The dependency on r would 
account for the vertical density and ve locity profiles in the ac- 
cretion column (see also iLeahvll2004alfbl) . We have verified that 
introducing a dependency of the emerging beam intensity on r 
does not have significant implications for the problem of model- 
ing the observed phase shifts. 

Our model of the radiation beaming comprises two emission 
components: one directed downwards and one upwards. Such 
beaming, at the end of our calculations, qua litatively reproduces 
the asymptotic beam patterns obtained by Sas aki et al.l {2011) 

7 The range of validity of this approximation was extended as de- 
scribed in Sect. [A] to account for the photons emitted above the NS 
surface in a downwards direction. 



with the decomposition method (see Eq. IA.91 l. The actual beam- 
ing is determined by the non-trivial interplay of special relativis- 
tic transformations due to the bulk motion of the emitting flow 
and anisotropic scattering. To estimate the magnitude of the spe- 
cial re lativistic effects, we adopt the column model of iBeckerl 
( 1998) for two values of their parameter e c (0.5 and 1) defined 
in their Eq (4.1). At the top of a 2 km column, the radiation ad- 
verted downwards is two to three times more intense than the 
radiation escaping upwards for an isotropic field in the matter 
rest frame. The escaped radiation becomes nearly isotropic at the 
bottom of the column, where the flow settles. According to this 
model, most of the radiation is emitted in the lower part of the 
column, below the radiative shock, which is predicted to reside 
at a few hundred meters above the NS surface. Weighting the 
downward and upward emission with their emissivity does not 
affect our main assumption, as we find that the total downward 
radiation is 40-70% more intense than the upward one. 

The bulk motion of the emitting plasma is not the only driver 
of the radiation beaming, because the scattering cross section 
are highly asymmetric in the strong magnetic field of HMXBs. 
Several stu dies have predicted a certain degree of anisotropy (see 
lKjausll200ll and references therein) for the radi ation emitted far 
from t he cyclotron resonances. More recently, Sch onherr et al.l 
(2007) have carried out MonteCarlo simulations and showed that 
the degree of anisotropy of the radiation in the energy range 
affected by the lines also depends on the optical depth of the 
plasma and does not exceed several percent (see their Fig. 7). No 
detailed studies have been carried out, to our knowledge, on the 
angular beaming of the radiation at energy close to the ones of 
the cyclotron scattering features, with the resolution required by 
the modern instruments. For this reason, our qualitative descrip- 
tion of the radiation beam through the introduction of upwards 
and downwards components was mainly inspired by the study 
of the rad iation spectral propert ies as a function of the viewing 
angle (see lSchonherr et al.l2007L a nd references therein). For the 
results of the existing modeling we notice that the absorption 
line is deeper in the direction perpendicular to the magnetic field, 
thus resulting in an effective suppression of photon number in 
that direction. These photons are scattered primarly along the 
direction of the magnetic field, where the cross section is sup- 
pressed, and would presumably lead to an enhancement of the 
beamed radiation along the column axis. As the Doppler boost- 
ing is not extremely severe in the velocity field of the infalling 
plasma, we conclude that a considerable fraction of the photons 
will be able to escape in the vertical direction and originate to an 
upward directed component. 

To summarize, we associate the downward component with 
the radia tion adverted by the plasma flow (as already proposed, 
see e.g.. lKrausll2001l) and the upward component, introduced 
here for the first time, with the (cyclotron-) scattered radiation. 
The latter component should be prominent close to the NS sur- 
face, where the flow is characterized by a reduced bulk motion, 
but gives rise to a non-negligible contribution also along the col- 
umn extension, owing to the relatively moderate Doppler boost- 
ing. We model the first component with a Gaussian centered on 
a = 150° and characterized by a width of 45° (normalization set 
arbitrarily to 100), and the second component with a Gaussian 
centered on the local zenith and characterized by a width of 45° 
and a variable normalization. Even though we use this heuristic 
approach, our model is able to qualitatively account for the most 



8 We concentrate our attention to the first harmonic, for which there 
is no photon spawning, which would add an unnecessary complication 
to our reasoning. 
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relevant properties of the intrinsic beam pattern emerging from 
the accretion column of 4U 01 15+63. 

Once Io,v(a) is given, the asymptotic beam pattern can be 
computed as a function of the angle between the axis of the col- 
umn and the line of sight to the observer (i^ b. s , see Appendix lAl 
for the details). In Fig. |9l we show the intrinsic beam patterns 
in the left hand panels, Io, v (a), for different values of the nor- 
malization of the upward component. The corresponding asymp- 
totic beam patterns are reported in the central panels. These were 
computed by assuming an accretion column of semi-aperture 4°, 
extended for 2 km above the surface of a canonic al NS with mass 
1 .4 M and radius 10km (column height from Nakaiim a et alJ 
2006). We notice that the shape of the beam is not very sensi- 
tive to the different parameters adopted for the column and the 
NS, provided that the upper part of the accretion column remains 
visible at large viewing angles (i// b$ ^ 150°). 

4.2. Simulated pulse profiles. 

Assuming the above asymptotic beam patterns, we computed the 
pulse profiles of 4U 0115+63 by using the geomet rical param- 
eters o f the system described in Figs. 2 and 3 of iRraus et alJ 
(1995, all measured with respect to the rotation axis of the NS). 
These are the direction of the line of sight to the source (©o), 
the colatitude of the two magnetic poles ©2), and their dif- 
ference in longitude (A). For the accretion column located on 
the /-magnetic pole, we can write the following relation between 
these angles: 

COS l/f b s = COS 0() COS @i + 

sin ©o sin ©,■ cos [Q(f - f ref ) - (i - 1)(jt - A)] , (2) 

where i = (1,2), and D. = 27r/P sp j n is the pulsar angular velocity 
measured at a reference time t =f re f Q We assume in t he following 
Qn = 60°, ®i = 148°, 2 = 74°, and A = 68° (ISasaki et all 
2011). As an example, we show the pulse profiles computed 
in the right hand panels of Fig. [9] by using the above geomet- 
rical parameters and the beam patterns reported in the central 
panels of the same figure. We did not consider other geomet- 
rical configurations, because we constructed the intrinsic beam 
pattern to reproduce the asymptotic beam pattern obtained for 
©o = 60°. However, we verified that phase lags can be modeled 
also for other allowed configurations by adopting a different in- 
trinsic beaming. 

By comparing Figs. Q] and [9] we note that the model de- 
veloped above can reproduce the shape of the pulse profiles 
of 4U 0115+63 relatively well at energies higher than the fun- 
damental cyclotron scattering feature. In our model, the bump 
in the profile that follows the main peak is due to the emis- 
sion from the accretion column closer to the observer's line of 
sight. However, we note that the double-peak structure of the 
upward emission (represented by the dot-dashed line at phase 
0.9-1.2 in Fig. [9j, which results from the model, does not seem 
to be present in the observed pulse profiles. The suppression of 
the flux between these two peaks (at phase ~ 1.1) is unavoid- 
able in the model, because it comes from the decreased flux of 
the asymptotic beam pattern at angles i/r b s < 30° produced by 
the self-obscuration of the column. The double peak could be 
avoided if the upward beam were emitted slightly outside the 
column (e.g., in a scattering halo near or on the NS surface), or 
if part of the radiation also escapes along the accretion column 

9 The spin period of accreting pulsars is not constant due to the torque 
exercised by the in-falling mass. 



axis in a pencil beam. Our simplified model could not reproduce 
this kind of emission, as this would have required either a ter- 
mination of the column at some height with emission from its 
top or the introduction of an emitting region around the column. 
However, we note that our calculation was not aimed at repro- 
ducing the exact shape of pulse profile, but rather at finding a 
plausible explanation for the phase lags observed in correspon- 
dence of the cyclotron lines, an effect which has never been re- 
ported before. 

Indeed, the phase lags measured from the data are also re- 
produced relatively well in the model (see Fig. [9]). According to 
our calculations, the displacement of the main peak in phase is 
due to the distortion of its shape at energies close to those of the 
cyclotron absorption lines, where the contribution of the upward 
Gaussian component to the intrinsic beam is enhanced. In our 
model, the larger estimated shift in phase is -0.05 (measured in 
phase units, for N2 = 75-100), and the corresponding value of 
the cross-correlation coefficient is 0.85. The analysis was limited 
here to the phases 0.0-0.5 that correspond to the mean peak of 
the pulse profiles, see right panels of Fig. [9] This is in fairly good 
agreement with the results found from the data. We also note that 
the location of the main pulse peak moves to lower phases at low 
energy (E < 5 keV), where a scattering halo on the NS surface 
dominates. This is a further indication that the scattered radiation 
is responsible for the observed phase shifts. 

Finally, we show that our model is compatible with the 
phase-resolved spectral properties of the absorption lines carried 
out in Sect. 13.31 In Fig. [TOl (upper panel), we report the viewing 
angles (^ bs,i) of the two accretion columns onto the NS with re- 
spect to the line of sight as a function of the phase. We also plot 
the average viewing angle (solid line) calculated for the com- 
bined radiation of the two columns (weighted for their relative 
contribution estimated by using the beam pattern in our model), 
which results in a relatively large angle (~ 120 - 150°) on the de- 
scending portion of the main peak (phases ~0.3-0.6), and signif- 
icantly lower values at other phases. To understand the meaning 
of these considerations, in the bottom panel of Fig. \TU\ we plot 
the relation between the angles a and iff, which are calculated 
from the relativistic ray-tracing of the model (see also Fig. lA.U . 
The former angle represents a good approximation of the rela- 
tive direction between the radiation emerging from the column 
and the local magnetic field, while iff is the trajectory deflection 
angle (see Fig. lA.U . which differs from the column viewing an- 
gle ((^obs) by at most a few degrees. We note that iff ~ 120 - 150° 
corresponds to a ~ 80 - 100°; therefore at the phases of the 
right flank of the main peak, the emerging radiation is emitted in 
a direction nearly perpendicular to the local magnetic field. At 
these inclinations, the resonant sc attering features are expected 
to be deeper and narrower (see ISchonherr et al.ll2007l and ref- 
erences therein), in agreement with the results found from the 
phase-resolved spectral analysis carried out in Sect. 13.31 

5. Discussion and conclusion 

In this paper, we analyzed archival BeppoSAX and RXTE ob- 
servations of 4U 0115+63 performed during the giant outburst 
of the source, which occurred in 1999. We investigated in partic- 
ular the shape of the pulse profiles of the source at energies close 
to those of the cyclotron absorption features (~ 12, 24, 36, and 
48 keV). At these energies, we measured a peculiar phase shift 
of the pulse profiles for the first time by using a cross-correlation 
method. 

To interpret these findings, we created a geometrical model 
for the emission from an accretion column of a neutron star 
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Figure 9. Left panels: intrinsic beam patterns (arbitrary units) 
emerging from the lateral walls of the NS accretion column as a 
function of the angle a with respect to the local radial direc- 
tion. The function is the sum of two Gaussians as described 
in Eq. (|A.9t . The value of the normalization of the upward 
Gaussian is indicated in each plot (A^)- Central panels: asymp- 
totic beam patterns (arbitrary units normalized to the values at 
170 degrees) as a function of the angle i/r b s between the axis 
of the accretion column and the line of sight to the observer. 
The beams are computed by assuming an accretion column with 
a semi-aperture of 4 degrees, and extended for 2 km above the 
surface of a neutron star with a mass of 1 .4 M Q and a radius of 
10 km. Right panels: synthetic pulse profiles computed for the 
geometry in 4U 0115+63 described in Sect. [4] The solid line 
corresponds to the total flux divided by two times its averaged 
value for plotting purposes. The dashed line represents the con- 
tribution from the farther pole from the line of sight to the ob- 
server, while the dot-dashed line indicates the contribution from 
the other pole. The gray vertical band corresponds to the loca- 
tion in phase of the pulse maximum. A leftward shift in phase is 
emerging while increasing the value of A^- 



that considers the most relevant relativistic effects. The intrinsic 
beam pattern was calculated by assuming a parameterization for 
the emission with both a downward and an upward beam, which 
mimics the presence of the radiation advected in the column by 
the in-falling accreting material and the effect of resonant scat- 
tering. This model permits the reconstruction of the asymptotic 
beam patterns from 4U 0115+63 as perceived by a distant ob- 
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Figure 10. Upper panel: inclination angles of the columns with 
respect to the line of sight to the observer obtained by assuming 
an intrinsic beam pattern with N2 = 25. The dashed line rep- 
resents the column that is farther away from the observer; the 
dot-dashed line, the other. The solid line is the average inclina- 
tion angle of the two columns weighted for their relative flux 
measured at the location of the observer. The dotted line is the 
synthetic pulse profile, scaled to fit into the plot. Lower panel: 
relation between the angles iff and a of Fig. lATI for different lo- 
cations from the center of the NS in units of the gravitational 
radius. The three locations correspond to the base of the accre- 
tion column, its midpoint, and the top of the column, respectively 
(calculated according to the geometry of the system defined in 
Sect. Hi. 



server and thus the calculation of the simulated pulse profiles 
from the source by adopting the geometrical parameters deter- 
mined bv lSasaki et"all(l201lh . 

Our method reproduced the general shape of the observed 
pulse profiles from 4U 01 15+63 relatively well and the proper- 
ties of their energy-dependent phase shifts described in Sect. 13. 21 
According to our interpretation, these phase shifts can stem from 
an energy-dependent beaming of the radiation, which we mod- 
eled as a variable contribution of the upward directed beam, 
whose origin can be ascribed to the resonant nature of the scatter- 
ing cross-section in a strong magnetic field. The modeled pulse 
profiles get distorted at the energies corresponding to the CRSFs, 
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creating a negative shift in phase that is in fairly good agreement 
with what is inferred from the observations. With the adopted 
system geometry, we were also able to naturally explain the 
phase-dependent measured properties of the absorption lines. 

Although we adopted a phenomenological description of the 
radiation emerging from the accretion column of an HMXB, our 
calculations revealed the key role of an energy-dependent beam- 
ing for the pulse profile formation mechanism. Based on this 
idea, we were able to point out a plausible explanation for the 
energy-dependent phase shifts that we discovered. A more com- 
plex approach, which will solve the radiation transport problem 
using the actual cross sections and a realistic column velocity 
profile, will provide a self consistent model of the spectral and 
timing characteristics of the XRB X-ray emission. 

Significant changes of the pulse profiles with the energy, 
especially close to the resonances of the cyclotron scattering 
cross-section, have been rep orted so far in a numbe r of differ- 
ent XRBPs (e.g., V 0332+53 fTsygankov et ai]|2006l) . Applying 
the analysis method developed here to these other sources is on- 
going and will be published elsewhere. 
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Appendix A: Calculation of the observed flux from 
an accretion column 

We present a relativistic ray-tracing computations of the beam 
pattern emitted from the accretion column surface of a slowly 
rotating neutron star to model the phase lags observed in the X- 
ray pulse profiles of 4U 01 15+63. We consider the geometrical 
properties of the column and the relativistic effects, i.e., light 
bending and the lensing effect in a Schwarzschild metric. Since 
the NS is slowly rotating, the relativistic travel time delay and the 
Doppler boosting have not been taken into account. The observer 
is located at positive infinity on the z-axis in the non-rotating ob- 
server frame with polar coordinates (r, 6, (p). The accretion col- 
umn reference frame has the polar coordinates (r',6',<f>') with 
the z'-axis coincident with the accretion column axis. See Fig. 
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Figure A. 1. The geometry of the systems (not in scale). The 
dashed line shows a photon trajectory from the accretion col- 
umn surface to the observer located at infinity, the dot-dashed 
line is the rotational axis of the NS. 




Figure A.2. Schematic representation of a photon trajectory 
(dashed line) from the accretion column surface to the distant 
observer for trajectories with a turning point (not in scale). The 
emission point is B and the angle to be computed is if/g. A is the 
point at which the trajectory reaches for the second time the ra- 
dial coordinate at which it originated, p indicates the location of 
the trajectory's periastron. 



lA.ll for the geometry of the systems. The deflection angle of a 
photon emitted from the accretion column surface is iff, while 
a indicates the angle between the local radial direction and the 
photon direction at the emission point. 

A. 1. Light bending 

The geodetic equation for an emit ting position of a p hoton de- 
pends on the impact parameter, b (Mis ner et alJ fl973). which is 
defined as 

b= J^^, (A.1) 
VI -rjr 

where r is the radial coordinate of the emitting surface element, 
and r s = 2GM/c 2 is defined as the Schwarzschild radius. 

The relation between a and the photon deflection angle, if/, 
i.e., the light bending, can be directly obtained using the approx- 
imated null geodet ic equation, strictly valid for a < n/2 (see 



1 a 

lBeloborodovl2002L for more details): 



1 - cos a - (I - cos 



*(.-*). 



(A.2) 



For each column surface element at a given column orientation, 
i^obs, (see Eq. 0, we determine the angle, a, and therefore, iff, 
at which the emitted radiation reaches the observer (the only tra- 
jectories considered for the pulse profile modelin g^) . 

For photons emitted towards the NS surface (a > n/2), the 
allowed a values can be estimated by imposing that the trajec- 
tory not be swallowed by the relativistic photosphere of the com- 
pact object. This translates into the relation tt/2 < a < a max , 
where 



= n - sin 









u 








r 


r t 



(A3) 



and r is the radial coordinate at which the photon is emitted from 
the column (e.g., point B in Fig. IA.2t . These trajectories have a 
"turning point" and thus possess a periastron; it is then conve- 
nient to express all relevant quantities in terms of the periastron 



10 We neglect values of the angle if/ such that if/ > n. For these val- 
ues multiple images of an emitting point could be formed. This is 
marginally relevant only for an accretion column located behind the 
NS, and thus not for 4U 01 15+63 (see Sect©. 



distance, p. For o- max > > tt/2, we first compute the im- 
pact parameter using Eq. ( IA.U . and then estimate the value p 
of the periastron by taking the largest real solution of the equa- 
tion p 3 - b 2 p + b 2 r s = 0. This equation is obtained by setting 



a = n/2 and r = p in Eq. (IA.U . as a p = n/2 at the periastron. 
Using Eq. (IA.2t . we find if/ p = cos~'(l - 1/ [1 - r s /p]) . The 
photon trajectory again reaches the radial coordinate r at point 
A (see Fip. IA.2t . where the angle formed by the photon with the 
local radial direction is &a = n — as- The angle fa can now be 
computed from Eq. iA.2\ . Finally, from simple symmetry con- 
siderations, 

fa = 2fa - fa . (A.4) 

This characterizes the considered trajectory fully, since the im- 
pact parameter for fa has the same value as for the deflection 
angle at the emission point, fa\. 

For each photon-emitting point we can now follow the pho- 
ton trajectory to verify if it is absorbed by the NS or accretion 
column surface. For a < n/2, the tra jectory can be dire ctly com- 
puted using the following equation ( Belo borodovl 2002 ) : 



r~ (1 - cos if/) b 2 
4(1 + cos if/) 2 sin 2 if/ 



r 2 (\ - cos if/) 
2(1 - cos if/) 



(A.5) 



where if/ values are in the range (0,if/). 

For a > n/2, we use Eq. (IA.U to link a to r* along the tra- 
jectory (the impact parameter b is constant in each geodesic). 
The corresponding value of if/ is found with the method outlined 
above. The trajectories that hit the NS or intersect one of the 
optically thick accretion columns are excluded from the compu- 
tation of the source flux measured at the observer's location. 

A.2. tensing and red shift 

We now compute the flux emitted by the column and measured 
by a distant observer. The flux emitted by a surface element from 
the accretion column, when it reaches the observed plane normal 
to the direction of the line of sight, is given by dF v = 7 v dQ, 
where the solid angle can be expressed in terms of the impact 
parameter (Mis ner et al H973h : 



dF v = 



I v (b, (f>)bdbd(f> 
D 2 ' 



(A.6) 
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Here, (p is the azimuthal angle around the line of sight of the 
plane containing the trajectory of the photon, and I v (b, <f>) is the 
differential intensity of the radiation emitted by the surface ele- 
ment when it reaches the observed plane at a distance D. 

The impact factor b depends on r and a, but not on <p. The 
angle a can be expressed as a function of r and if/, as shown 
in the previous section. We further note that the differential in 
Eq. ( IA.6t can be expressed as db = (db/dr)^dr + {db/dif/) r dif/, 
and is computed along a line on a geometrical cone at constant 
<f>, which is by definition a radial line. Therefore we can simplify 
the differential and write db = (dbfdr)^,dr, since if/ is constant 
along a radial line. 

Up to this point we expressed the equations in the observer's 
reference frame defined above. To integrate the emission on the 
column surface, we transform the integration variables to the col- 
umn's frame to use quantities directly related to the accretion 
column. The Jacobian of the coordinate transformation involves 
only the variables r (constant) and <p, and is: \d(f>/d(p'\; therefore, 
Eq. dA.6t can be written as 



dF v = 



I v {a,<p',r)b{a,r)f r \^\d<t>'dr 
D2 ' 



(A.7) 



We still need to express I v in the system of reference of the col- 
umn to account for the gravitational red-shift. Since 7/v 3 is a 
Lorentz invariant, we can write 

Iy = (jf k,v = (l - /0,v - (A.8) 

As already described in Sect.|4] knowing the exact functional 
shape of Io, Vo (f, a) would require a detailed treatment of the cy- 
clotron scattering process in a strong magnetic field (B ~ 10 12 G) 
and of the relativistic beaming caused by the bulk motion of the 
plasma in the accretion column above the surface of the NS. 
This is outside the scope of the present paper. Instead, we aim 
at reproducing the beam pattern in the energy range of interest 
(> lOkeV ) found by u s ing th e pulse decomposition method de- 
scribed in Sasa ki et alj For this purpose, we use a com- 
bination of two Gaussians, one pointed downwards and one up- 
wards characterized by different amplitudes: 




Figure A.3. Column beam pattern seen by a distant observer and 
calculated for an isotropic emission in the case of a NS with a 
radius r - 3rs. The accretion column is assumed to have an half- 
aperture of 5 degrees, and an heights of (from bottom to top) 
3.03, 3.15, 3.3, 3.6, 3.9, 4.2, and 4.6 rs- The angle if/ is measured 
from the line of sight to the axis of the column. 



column for its total obscuration in the antipodal position should 
be 3.94r s . With the approximation used in our calculations, we 
obtained a value of 3.85r s . This discrepancy can only be rele- 
vant for a height of the column that is close to this critical value 
and for angles close to the antipodal position. Our approxima- 
tion thus gives a satisfiable result for the accuracy required in 
the present work. 



I (a, r) = 2_j N > ex P 





la-ai\ 




\ <ri)\ 



(A.9) 



In the equation above, a x = 150°, a 2 = 0°, cr x = 45°, <x 2 = 45°, 
N\ = 100, and No assumed values from to 100 to mimic a 
variable contribution from the cyclotron scattered radiation (See 
also Sect.g}. 

The total flux emitted by the column is finally estimated from 
Eqs. (IA.7t . 1A.81 . and ( IA.91 > performing a Monte-Carlo integra- 
tion with <p' e (0, 2n), and r e (r^s, '"max)- The maximum height 
for the accretion column is equal to r max - r^s • We did not in- 
clude in the total flux the contribution of the photons that hit the 
NS surface and/or are absorbed by the column. 

To verify the consistency of our results with similar calcula- 
tions published before, we show in Fig. IA.3l the beam patterns 
obtained from an isotropic emission on the lateral surface of an 
ac cretion column with the sa me configuration adopted in Fig. 3 
of iRiffert & Meszarosl ([988). We found fairly good agreement 
between our results and those published by these authors. There 
is only a negligible discrepancy in the estimated fluxes of an ac- 
cretion column extending up to 3.9 gravitational radii above the 
NS surface. Exact calculations predict that the critical height of a 
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